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Many nonlinear PDE display a coarsening dynamics, i.e. an emerging pattern whose typical length 
scale L increases with time. The so-called coarsening exponent n characterizes the time dependence 
of the scale of the pattern, L(t) « t n , and coarsening dynamics can be described by a diffusion 
equation for the phase of the pattern. By means of a multiscale analysis we are able to find the 
analytical expression of such diffusion equation. Here, we propose a recipe to implement numerically 
the determination of D(\), the phase diffusion coefficient, as a function of the wavelength A of the 
base steady state uo(x). D carries all informations about coarsening dynamics and, through the 
relation \D(L)\ ~ L 2 /t, it allows to determine the coarsening exponent. The main conceptual 
message is that the coarsening exponent is determined without solving a time-dependent equation, 
but only by inspecting the periodic steady-state solutions. This provides a much faster strategy 
than a forward time-dependent calculation. We discuss our method for several different PDE, both 
conserved and not conserved. 



I. INTRODUCTION 

In this paper we will focus on the nonlinear dynamics 
of one dimensional Partial Differential Equations (PDE) , 
having the form [l| 



d t u = Af[i 



(1) 



where Af is a general nonlinear operator, not depending 
explicitly on time t and space x and whose trivial solution 
u = has a linear stability spectrum whose modes of suf- 
ficiently small wavevector are unstable. More precisely, 
setting u(x,t) = uexp(iqx + uit), and keeping only terms 
which are linear in u, one obtains a dispersion relation 
relating 10 to q. Two well-known dispersion relations [lj 
are, 

w{q) = 1 — q 2 , nonconserved models (2a) 

ui(q) = q 2 — q 1 ■ conserved models (2b) 

In both cases one has that u}(q) > for q < 1. In 
other words, the trivial solution u = is unstable for 
q < 1. This means that modes with arbitrary small 
wavenumbers are unstable, and therefore, in principle, 
spatial structures with large wavelength can develop, so 
that the notion of coarsening makes sense. This implies 
that we exclude situations where where there is a small 
q cut-off, that is, when there is a minimal wavenum- 
ber below which the trivial solution is stable. A typical 
example in this category is the Swift-Hohenberg equa- 
tion PJ , for which the dispersion relation takes the form 
uj(q) = a+(q — l) 2 , with a a real parameter. This disper- 
sion relation means that if a < 0, then there is a minimal 



wavevector q m i n = 1 — y/a below which the trivial solu- 
tion is stable. We do not expect in this case a structure 
with a wavenumbcr smaller than q m in to take place (no 
perpetual coarsening is a priori expected). While this 
statement complies with intuition, we must keep in mind 
that the stability evoked above follows from a linear anal- 
ysis. Therefore, since the equation of interest are nonlin- 
ear, exceptions are not excluded 0. Having in mind the 
occurence of some exceptions, we will focus on systems 
having a linear dispersion relation of the form (|2|) . This 
holds for many nonlinear PDE. The most known equa- 
tions are the Ginzburg-Landau (GL), the Cahn-Hilliard 
(CH), the Kuramoto-Sivashinsky (KS) and their variants 
or generalizations [l[ . Equations of this class will be in- 
troduced and discussed in the next Section. 

Note that having unstable modes for small q as in @ 
is not a sufficient condition. A prominent example is the 
KS equation, which shows spatio-temporal chaos, with 
wavelengths of constant average size being continuously 
destroyed and created. This is not surprising: the anal- 
ysis of the linear stability spectrum does not univocally 
determines the nonlinear dynamics. 

One important feature that we will discuss here is the 
determination of the branch of periodic steady-state solu- 
tions, which acquire a special status in our investigation 
of coarsening dynamics. This is a priori not obvious, 
since coarsening is a time-dependent process, but it has 
been shown before that stationary configurations of pe- 
riod A, uq(x), play a major role in determining the type 
of nonlinear dynamics: for a certain class of nonlinear 
equations in ID, two of us have proven Q that coars- 
ening is possible if and only if the wavelength A is an 
increasing function of the amplitude A of the solution 
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uq(x), X(A) > 0. The presence of a maximum in the 
curve X(A) signals the phenomenon of interrupted coars- 
ening, while a decreasing X(A) signals no coarsening at 
all. Finally, for some equations like KS, \(A) is not a 
single-value function and the curve displays a so-called 
turning point. 

All these examples show that steady states contain im- 
portant pieces of information. However, how is it possi- 
ble that dynamics depends on stationary properties only? 
The answer to this question is contained in the concept 
of phase instability. For example, in ID, a steady-state 
periodic solution uq(x) has a wavelength (or wavenum- 
ber q = 2ir/\) A such that uq(x + A) = uq(x). In this 
case the phase cj> of the pattern is </> = qx. However, a 
perturbation of the wavelength of the pattern (or equiv- 
alcntly a perturbation of the phase </>) may reveal insta- 
bility (called phase-instability), meaning that the wave- 
length will evolve in time, and if this is the case for any 
q, one may expect coarsening. A global shift of a peri- 
odic pattern in the x direction by a constant value xq 
is a also a solution, i.e. uq(x + xq) = uq(x), owing to 
translational invariance (also called a Goldstone mode). 
This constant shift corresponds to a phase perturbation 
of infinite wavelength, and since it is a stationary solu- 
tion, it has an infinite relaxation time towards the orig- 
inal solution uq(x). In other words, a large wevelcngth 
perturbation of the pattern is expected to have a large 
time scale. This entails that the most dangerous pertur- 
bations are those of longwavelengths of the phase of the 
pattern. To this notion explicit, we introduce a small 
parameter e expressing the fact that the phase modes 
of interest have a small spatial gradient, on the scale of 
the wavelength of the pattern. An appropriate way to 
deal with this question is to introduce a multiscale anal- 
ysis: the fast spatial scale x (where the periodic profile 
of pattern varies on a scala of order unity) and a slow 
scale X (expressing the longwavelength perturbation of 
the phase of the pattern). The slow spatial scale, implies 
also a slow time scale (which can also be inferred from 
the above discussion about large relaxation times) T. It 
turns out that (see later) X = ex and T = e 2 t. 

If <p designates the actual phase of a pattern (not nec- 
essarily a periodic one, but some perturbed pattern), it 
will depend on space and time, e.g. because Uq{x) is per- 
turbed, and acquires a slow spatial and time dependence 
(i.e. it depends on X and T), which can be described by 
a phase diffusion equation (see next Section), having the 
form 



d T ip = D(q)d X xip, 



(3) 



where ip,X,T arc suitably rcscalcd versions of (f>,x,t, 
and the phase diffusion coefficient D(q) only depends 
on properties of the steady state solution of wavelength 
A = 2n/q. A negative D(q) signals a phase instability, i.e. 
coarsening. Other scenarios are possible Q. Even if for 
several nonlinear equations it is feasible to obtain the an- 
alytical form of D(q) in the limit of large A, in the generic 
case it is not possible and we must determine it numeri- 



cally, by computing the steady-state solution uo(x), from 
which we can state if the pattern is stable (no coarsen- 
ing), or unstable (coarsening). We will then show how 
the coarsening exponent (in the temporal power law) can 
be extracted from these considerations (i.e. without any 
time integration of the equations). 

In this paper wc start shortly describing the multi- 
scale method and apply it to several ID PDEs, there- 
fore extracting the quantities we need to determine D(q) 
numerically. We will see that we need two functions: 
(i) uq(x), the stationary solution, which satisfies the 
equation Af[uo] = 0, see Eq. ([T]). (ii) Vq(x), solution 
of the equation C^[vq(x)} = 0, where C is the Frcchct 
derivative of Af and is its adjoint operator. We will 
show how it is possible to implement the numerical deter- 
mination of such functions, therefore the determination 
of D(q) and, finally, how we can get the coarsening ex- 
ponent n, describing the time dependence of the typical 
scale, L w t n , when coarsening is present. Particular 
attention will be devoted to the conserved Kuramoto- 
Sivashinsky (cKS) equation, because in such case C is 
not self-adjoint and special care must be devoted to the 
determination of vq(x). 



II. FROM THE ID NONLINEAR EQUATION 
TO THE PHASE DIFFUSION EQUATION 

In this Section wc introduce all nonlinear equations 
we are going to discuss, we give a very brief sketch of 
the multiscale approach used to extract the phase diffu- 
sion equation, and we finally provide the expression of 
the phase diffusion coefficient D. In some cases, we just 
report existing results, but in the case of the cKS equa- 
tion, we provide the original derivation of D, which will 
be seen to deserve special attention. 



A. The models 

Let us start by listing all equations of interest, the 
nonconserved first: 

dtu = u — u 3 + d%u, GL equation (4a) 

u 



d t u 



cLu, a- models 



(1 + u z ) 

then their conservative counterparts: 

d t u = —dl [u — u 3 + d x u] , CH equation 



(4b) 



d,u = -at 



» r,2 

(1 + v?) a xl 



, ca-models 



(5a) 



(5b) 
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and, finally, the conserved Kuramoto-Sivashinsky equa- 
tion 

d t u = —d 2 [u — rd x u + d 2 u — (d x u) ]. cKS equation 

(6) 

GL and CH are classical equations @ describing, e.g., 
phase separation processes for a scalar order parameter, 
either conserved (CH) or nonconserved (GL). The ca- 
model for a = 1 firstly appeared in problems of crys- 
tal growth Q and starting from that specific conserved 
model it has been natural to extend it [8| to a ^ 1 and 
to the nonconserved case. Finally, the cKS equation ap- 
peared in the study of sand-ripple dynamics |9j , then in 
step- bunching phenomena • For r = it also appears 
in the study of wandering cristalline steps [ill, EH ■ 

Eqs. (IH have the general form dtu = A{u) + d 2 u and 
their linear dispersion relation has the form (f2a|) . All 
conservative equations, (0 and ([6]), have the form dtu = 
— d% (■ ■ ■ ) and their linear dispersion relation is like (|2bj) . 
For Eq. (jSJ), that is true for r = only. However, the 
r-term can be cancelled out with a tilt transformation: 
u — > u . + tx/2. We keep it here for sake of generality. 

All previous models, conserved and nonconserved, dis- 
play perpetual coarsening, because, as analyzed in the 
next subsection, they have a branch uq(x) of steady 
states whose X(A) characteristic is an increasing, diverg- 
ing function. 

B. The multiscale method and the phase diffusion 
equation 



This is enough, along with the e-expansion of the so- 
lution, u = u + eui, to perform a correct perturbative 
analysis of Eq. (p}. Using Eqs. (JTHHJ) , the general non- 
linear equation (TTJ) rewrites as 

e(dd,uo)d T ip =(A/"o + eA/l)[wo + 

(9) 

=M [u ] + eC [ui] + eA/i[uo], 

where Cq is the functional, Frechet derivative of AV 
Since Mq[uq] trivially vanishes, we are left with the first 
order terms which can be rearranged, so that the function 
uq appears in the source term, g, of the linear differential 
equation for ui, 

C [ui] = (d^uo)d T il} —Mi[uq] = g. (10) 

According to the Fredholm alternative theorem [l4|; g 
must be orthogonal to the kernel of the adjoint operator, 
Cq. In conclusion, if 

4w=o, (ii) 

then 

(«0,(^uo)Or^-M[«o])=0. (12) 

In the next subsection we are going to see that Eq. ([12")) 
can be recast in the form of a phase diffusion equation, 

d T i> - D(q) d X xi>, (13) 



We are now going to discuss the multiscale method, 
with a focus on the cKS eq. This step has two motiva- 
tions. First, cKS eq. is a completely non trivial equation, 
whose nonlinear dynamics has been studied only through 
direct numerical simulations or through a simplified par- 
ticle model. We are not aware of analytical nonlinear 
studies. Second, its treatment with the multiscale ap- 
proach is representative of the method, when the Frechet 
derivative is not self-adjoint. 

The basic idea of the multiple scale expansion 
method [l3| is to introduce slow space and time variables, 

X = ex, T = e 2 t, (7) 

along with the fast space variable, x. No fast time vari- 
able exists, since the base state uq(x) is a stationary 
state. Identifying the fast phase <j)(x, t) = qx and the 
slow phase xp(X,T) = e<fi(x,t), where q = q(X,T), we 
can rewrite the space and time derivates, 



where the phase diffusion coefficient, D(q) is a nontrivial 
functional of uq and vq. 



C. The phase diffusion coefficient 

The application of the multiscale method to Eqs. (|4| 
and (|5|) gives, respectively 

d g (q (O^uq) 2 ) 

D(q) = — j r— — , nonconserved Eqs. (jlj) (14) 

((<Vo) 2 ) 

and 

D(q) = q 2 x 2 '-. conserved Eqs. |S) (15) 

\ u o) 



d x = qdtf, + ed x = qdj, + eipxxdq, (8a) 
d t = eV-TcV (8b) 



For the cKS equation ©, we give here below a few 
details for the derivation of D(q). In that case, the oper- 
ators appearing in Eqs. (jWT^j) are 
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FIG. 1: (a) Evolution of the residue poo of the Catm-Hilliard equation <|5a|) as function of time for different values of At (reported 
in the legend). The initial condition is a sinusoid discretized by N = 512 points with A = 12.8, corresponding to q = 0.49087. 
Note that for a fixed At the residue saturates at an asymptotic value which depends on the magnitude of the time steps, (b) 
Evolution of p x as function of the iteration step m for the same simulation condition of the left figure. The black points are 
the value of poo for a numerical integration of Eq. (|5a[) with a variable time stepping. Starting from At = 0.025, the magnitude 
of At is reduced of a factor ten (St = 0.1) every time the condition Apoo < 10~ 4 is verified. The numerical integration stops for 
Poo < 10~ 9 , value reached for At = 2.5 x 10~ 6 . Note that the sudden improvements of the accuracy of the estimated solution 
tto are due to the reduction of At. 



•M) [«o] = -Q 2 dj [(1 - rqd^ + q 2 dl) u Q - q 2 (d^Uof 
Co = -q 2 d 2 [1 - rqd* + q 2 2 - 2q 2 (d^ ) fy] Ui, 
M [uo] = -ipxx {q 2 dl [(2qd q + 1) d<j,UQ - Td q u - 2q (d^uo) d q u ]} , 
4[v] = -q 2 {l + rqd^ + q 2 d 2 + 2q 2 [(d 2 u n ) + (9^ ) } d 2 v. 
Finally, from Eq. (|12[) we can determine the diffusion coefficient, 

D (l) = -Q 2 (dlv , (2qd q + 1)8^0 - vd q u - 2q(d q u ) d<f,u }/(vo, d<j,u }. 



(16) 
(17) 
(18) 
(19) 

(20) 



When Co is self-adjoint, as for Eqs. (j4]), traslational 
invariancc immediately provides the result vq = d^UQ. 
For Eqs. ©, £j ^ Co, but Eq. 1^) can be easily solved 
analytically and vo is still related to uo through some dif- 
ferential operator. In all these cases, covering Eqs. (|4ajl - 
(|5b[) . D can be expressed as a function of uq and its 
derivatives only, see Eqs. (fT4")) - (fTl)j) . Finally, there are 
cases where C\ ^ Co and vo cannot be determined ana- 
lytically. The cKS equation is a typical example. In the 
next Section we discuss the most general implementation 
for a numerical determination of uo,vq and D. 



III. NUMERICAL IMPLEMENTATION 



The numerical procedure we have developed to find 
the value of D for a given A can be split in three parts: 
(i) The determination of the stationary solution uq by 
employing a pseudospectral algorithm and, only if re- 
quired, a Newton-Raphson method to increase the accu- 
racy of the output of our pseudospectral code, (ii) In the 
second step, the function vq is found by discretizing the 
adjoint £j of the Frechet derivative of 7V"o with high or- 
der finite differences and finding the kernel of the sparse 
matrix M containing the coefficients of the discrete ver- 
sion of Cq. Clearly, we skip this step when we are able 
to obtain vo analytically, (iii) Lastly, the value of D(\) 
is computed with a standard numerical integration after 
the evaluation of the derivatives (in <j) and q space) of uq 
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%I2 Jt 3n/2 2jc tc/2 Jt 3n/2 2% 

§ 4> 

FIG. 2: (a) Normalized stationary solutions uo of equation (|23|) with r = for different values of q. The differential operators 
are discretized in an uniform gird with A<f> — 0.025. For visualization purposes we have chosen a normalization constant 
A = max |it(0)|. (b) Solutions of the adjoint problem £q[uo] = 0, see equation (|19[) . with r = 0. As for the left figure, the 
spatial discretization is A<j> = 0.025. 



and, if necessary, v (only for the cKS Eq.). These three 
steps are iterated for each A and, finally, the coarsening 
exponent is estimated using the relation \D\ ~ L 2 /t. 



the numerical implementation). For a comparison be- 
tween different pseudospectral algorithms see [2(| and 
references therein. 



A. Determination of the stationary solutions 

The time dependent problem (|TJ) can be integrated very 
efficiently by using a pseudospectral algorithm. These 
methods are widely used for the integration of nonlin- 
ear PDEs because they combine the estimation of spa- 
tial derivatives in Fourier space with the calculation of 
nonlinear terms in real space, hence avoiding the compu- 
tational overhead of the convolution step [15l - tl8j . Pseu- 
dospectral methods mainly work in Fourier space and 
their strength resides in their high accuracy for smooth 
solutions. 

It is convenient to restate the nonlinear problem (TTJ) 
in <fi <G [0, 27r] space, by replacing the derivatives d x with 
qdfi. In this way, stationary solutions obtained at differ- 
ent values of q are mapped into the same interval so that 
they can be directly compared. 

Every nonlinear PDE such as Eq. ([T]) can be divided 
into a linear operator ui, the dispersion relation, and a 
nonlinear operator N, such that in Fourier space Eq. (TTJ) 
reads 



d t u k (t) =w k u k {t)+N[u{t)} 



A- ' 



(21) 



where Nfit]*, is the Fourier transform of the nonlinear 
part of the PDE, and the discrete wavevectors k belong 
to the interval [-N/2,N/2], with the same Ak = 1 for 
every A (here N is the number of points used to dis- 
cretize the <fi interval). Eq. (1211) i s evaluated by using 
the integrating factor technique [I5l.[l8j and time stepped 
with a fourth order Adams-Bashforth-Moulton predictor- 
corrector scheme (IFABM4) (see for the details of 



The IFABM4 method has been used in this paper to 
study the time evolution of Eq. ([I]) and to characterize 
the coarsening exponent by measuring the mean wave- 
length of the emerging pattern, L(t), as a function of 
time t. In order to obtain a reliable estimate of n, the 
function L(t) must be averaged over several realizations, 
each one with a different initial condition. 



Stationary solutions of nonlinear PDEs are normally 
found with specialized approaches that have better per- 
formances compared to any numerical integration of the 
dynamics of Eq. ffl. Although Point Collocation or 
Multigrid methods [19| are less computationally expen- 
sive than our IFABM4 code, we decided to rely on the 
same algorithm to highlight the advantages of a strat- 
egy based on the the phase-diffusion equation. More- 
over, even if it is well established that the integrating 
factor technique can lead to a wrong estimation of the 
fixed points of nonlinear ODEs [2l| and more accurate 
algorithms have been already developed for ODEs (2~l| 
and PDEs [2(| like Eq. (TTJ), a small modification of the 
IFABM4 algorithm allowed us to find the stationary so- 
lutions uq with a reasonable precision (enough for our 



purposes). 

Starting from an initial guess function u(x, 0), the dy- 
namics of Eq. (|2"Tj) leads to the closest stationary solution. 
Obviously this can happen only when al least one uq ex- 
ists for the fixed value A and for the chosen parameters 
(a, r, . . . ) considered during the integration of Eq. (|21l) . 
The error committed in the estimation of the stationary 
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FIG. 3: Left vertical axis: Stationary solution uo (thick full line) of Eq. (|4b[l with a = 2 and A = 102.4, corresponding 
to q = 0.061359. The parameters of the IFABM4 algorithm are N = 4096, At = 0.1, St = 0.1, Ap min = 10~ 4 , and a 
stopping condition poo < 5 x 10 -9 . Right vertical axis: residue p(cf>) before (blue, thin line) and after (red, dashed line) the 
Newton-Raphson minimization. 



solution is easily quantified by inspecting the residues 

P {4>) = r'w^ + N^)], 

Poo = max \p{4>)\, (22) 

during the numerical integration of the PDE. As shown in 
Fig. for a fixed value of At we observe that p^ satu- 
rates at an asymptotic value depending on the magnitude 
of the time steps. Even if starting the numerical integra- 
tion of Eq. (|21"j) with a tiny At ensures a small asymptotic 
residue, this will increases enormously the computational 
cost required to find ito, therefore impairing the applica- 
bility of this method to many practical cases. 

A possibility to overcome this trade-off between accu- 
racy and computational effort arises when we take into 
account that here our goal is not to follow the real dy- 
namics of Eq. ([21"]) but rather to obtain a reliable estimate 
of its stationary solutions. The numerical convergence 
of u(t) to uq can be increased by reducing dynamically 
the value of At. Actually, by monitoring the residue 
and its variation Ap™ = \p™ - p^r m "\l ' pZ (here the 
superscript m stands for the discrete time) every prede- 
fined number of iteration steps m c , we are able to sense 
when the dynamics of the PDE has reached the station- 
ary state for a given At. In fact, a drop of Ap™ below 
a given threshold Ap rn i„ means that the residual error 
cannot be diminished anymore, signaling that is the ap- 
propriate moment to reduce At by multiplying its value 
for a factor 6t < 1. After every reduction of the time 
step, the predictor-corrector scheme has to be warmed- 
up by a method of equal accuracy, such as a fourth order 
Runge-Kutta (see |X] for its numerical implementation) . 



In general, the values of Ap m i n and St have to be inferred 
empirically for each problem: for our analysis, we have 
used Ap m i n = 10 -4 and 6t = 0.1. Additionally, m c must 
be a large number in order to limit the computational 
overhead required for the estimation of the residue and 
its variation. For one-dimensional problems a value of 
m c between 10 3 and 10 4 is a reasonable choice. As an 
example, in Fig. we report a comparison between the 
constant and the variable time stepping in case of the 
Cahn-Hilliard equation (|5aj) . 

In some cases the solution found by the pseudospec- 
tral algorithm is not accurate enough to ensure a good 
estimation of the coarsening law, through the evaluation 
of D{q). This usually happens for large system sizes, so 
that the algorithm converges slowly, or when the vari- 
ations of the stationary solution are extremely localized 
and the resolution of the numerical method is insufficient 
to resolve these jumps (spectral methods are not suit- 
able to handle discontinuities like shocks). To improve 
the residue of u we have employed a Newton-Raphson 
(NR) method applied to the finite difference discretiza- 
tion (in <j> space) of the equation M[u] = 0. To be more 
precise, in case of an operator M that can be written 
as a second order derivative of another nonlinear oper- 
ator M r (that is M[u] = — d 2 M r [u}) we have found the 
solutions of equation A^it] = and then subtracted to 
these functions their mean value (uq). In fact, these type 
of nonlinear operators are conservative so that the mean 
value of u is preserved by the dynamics of the equation, 
and, additionally, we can deal with an operator that is 
less stiff than the original one. Besides, it is sufficient to 
find the stationary solution of the GL equation (f4"a|) to 
compute the phase diffusion coefficient of the CH equa- 
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FIG. 4: (a) Phase diffusion coefficient D estimated with reduced expressions (|f 4[) and (|f 5[) for the Ginzburg-Landau equation 
(|4a|l (blue circles) and for the Cahn-Hilliard equation (|5a|l (black rhombi), respectively. Solid lines represent values of D 
calculated from Eqs. (|33[) and (|34|) . Note that for these two cases D is negative until the sensibility of the algorithm (in this 
case of order 10 -7 ). The stationary solutions uq have been obtained with a mesh size Arc = 0.025, a varying time step At 
decreasing from 1 to 10 -6 , and a maximum residue p m = 10 -8 . The q derivatives in (|14l) and (|15l) are estimated by evaluating 
four equally spaced wo around q = 2-k/\ with a Aq = 2.5 x 10~ 3 . (b) Estimation of the coefficient C [A = C" 1 ln(t)] for the 
Ginzburg-Landau equation (|4ap (blue circles) and the Cahn-Hilliard equation (|5a|l (black rhombi). This coefficient is estimated 
through a linear regression of In (|D|/A 2 ) by applying a sliding window binning four consecutive values of \D\ from the left 
figure. The wavelength A is the average of the different wavelengths of the points binned together. All units are arbitrary. 



tion (|5a|) (the same for the a-models). In the following, 
as an example, we detail how we have implemented the 
Newton- Raphson method to find the stationary solutions 
of the conserved Kuramoto-Sivashinsky equation. 

The dynamic evolution of the cKS equation presents 
several numerical difficulties, which can be ascribed to 
the nature of its stationary solutions: a sequence of arcs 
of parabola whose amplitude grows quadratically with A, 
while their joining regions have a diverging curvature and 
a vanishing size [9|, Il2| . This leads to the blow up of the 
integration scheme for large time steps. More specifically, 
the numerical integration of a system with Ax = 0.1 and 
N = 2048, i.e. L = 204.8, must be performed with our 
IFABM4 code by choosing At — 10~ 5 in order to prevent 
the blow up at t max ~ 10 3 . For this reason, the phase 
diffusion method is very suitable to estimate numerically 
the coarsening law for the cKS equation. In this case, 
instead of performing a first estimation of the stationary 
solutions uq by means of the pscudospcctral algorithm, 
we have solved directly the equation 



where A4 , Bi are given by 



(1 - rqd^ + q 2 dl) u - q 2 {d^uf = 0, 



(23) 



by discretizing the differential operators with six order 
centered finite differences. In this way Eq. (|23j) is trans- 
formed in a set of N nonlinear equations in the space 
of the variables Ui (with i = 1, . . . , N) that is solved by 
finding the zeros of the set of functions 



-4, 



1 3 



j'=i 



Bi = (d% 



A0 2 



boUi + ^2 bj {u i+] + m-j) 



(25) 
, (26) 



and the coefficients aj,bj are listed in Tab. |TJ These vec- 
tors are built by taking into account the periodic bound- 
ary conditions of Uj. The Jacobian J of this system of 
equation is needed in order to solve the problem. The 
elements of this matrix are Jij = d Uj Fi, which is a band 



matrix with 



J, 



F, 



q 2 B, - qAi (r + qAJ 



(24) 



(A</>) 2 ' (27 ) 
J ht±J = (A0)" 2 [q 2 b 3 T qaj (rA0 + 2qA t )] , 

where j = 1,2,3. Through the Matlab® func- 
tion fsolve, the values of Ui for which Fi = 
are readily found by means of the standard large- 
scale method implemented in the package, i.e. the 
trust-region-ref lective algorithm (for details sec 
[22I [23( 1 and the Matlab® documentation). Some sta- 
tionary solutions Mo of the cKS equation for t ~ are 
reported in Fig. [5k.. 

To conclude this section we present one case in which 
the Newton-Raphson method has been used to decrease 
the residue of the stationary solution found by our 
IFABM4 code. As working example we consider the non- 
conserved a- model (|4bp with a = 2. In Fig. [3] we show 
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FIG. 5: (a) Modulus of the phase diffusion coefficient of the non-conserved a-models (|4b[) for different values of the parameter 
a. The stationary solutions have been obtained through two steps: i) a first approximation of uo is obtained by employing 
the IFABM4 algorithm up to a precision of about p m ~ ICF 7 , ii) the final value of the stationary solution is found by refining 
the result of the previous step through the Newton-Raphson method. We used a spatial discretization A(f> = 0.025 and 
Aq = 1.0 x 10~ 4 for a = 1, 2, whereas A0 = 0.05 and two different Aq = 1 .0 x 10~ 4 , 2.5 x fO" 5 ensuing a ratio Aq/q < I0" 2 
for the case a = 3 (see right figure). The coefficient D has been estimated through equation (114[) . According previous 
analytical calculations [ij its asymptotic value saturates to a constant value, i.e. \D\ ~ 1, when a < 2. The solid line is a fit 
\D\ = ag/(a\ + \n(q), showing the logarithmic correction for a — 2. (b) Estimation of \D\ by means of (| 15j) for the conserved 
a-models ()5b|) with a — 1,2,3. We used the same stationary solutions found for the non-conserved a-models. The expected 
power law for these conserved models is \D\ ~ q 2 [J]. The blue solid lines are guides to eyes with slope equal to two. All units 
are arbitrary. 



the solution uq (black line) and the residue p(4>) before 
(blue, thin line) and after (red, dashed line) the NR 
minimization. The residues of the starting guess (blue 
line) had Poo = 2.14 x 10~ 9 and p 2 = 1.35 x 10" 10 , but, 
after the minimization step, these residues reduced to 
Poo = 1.34 x 10" 9 and p 2 = 1.12 x 10" 11 . Note the dras- 
tic reduction of p(cf>) in the regions characterized by small 
values of d<pUo leading to a p 2 that is one order of magni- 
tude smaller than before the Newton-Raphson minimiza- 
tion. 



index i 





1 


2 


3 


4 


di 




3/4 


-3/20 


1/60 




bi 


-49/18 


3/2 


-3/20 


1/90 




a 




-61/30 


169/120 


-3/10 


7/240 


di 


91/8 


-122/15 


169/60 


-2/5 


7/240 



TABLE I: Coefficients of the six order finite difference dis- 
cretization of the differential operators used in the Newton- 
Raphson minimization. 



B. Solution of the adjoint operator and calculation 

of D 



account that in Eq. ([19]) we have derivatives up to fourth 
order, first and second derivatives arc discretized accord- 
ing to Eqs. (|25l) . (f2"6]l . and third and forth derivatives are 
given by 



(3 



(3 



(A0)- 



(A0)" 



4 



(v 



i+j 



d Q Vi 



J'=l 



(28) 



(29) 



where the coefficients Cj and dj arc listed in Tab. U] 

The derivatives 9?uo and O^uq are computed in Fourier 
space. The elements of these vectors multiply the con- 
stant coefficient of the Finite Difference discretization of 
the differential operators of Eq. (fl~9|) . thus forming the 
sparse matrix M, which represent the discrete version of 
Eq. (fTrj|) . Then, by means of the Matlab® function eigs 
(that is based on the Fortran library ARPACK), we search 
for the eigenvalues of M which are in modulus closer to 
zero and their associated eigenvectors. The eigenvector v 
related to the smallest eigenvalue is finally normalized in 
order to have a solution vq with zero mean and amplitude 
equal to one 



The solution of the adjoint problem Cq[v] = is easily 
found through the finite difference discretization intro- 
duced in the previous section. Again, let us consider the 
conserved Kuramoto-Sivashinsky equation. Taking into 



t'o 



max(|{)|) 



Some functions v for the cKS equation for r 
reported in Fig. 



(30) 
are 
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FIG. 6: (a) Estimation of the power law behavior of the absolute value of the phase diffusion coefficient of the conserved 
a-models for different values of a. The exponent 7 (that is \D\ ~ g 7 ) is computed by a linear regression of three consecutive 
points of each curve of Fig. [5b. The center of the sliding window q is the mean value of the abscissa of these points, (b) Relative 
error of the estimated exponent 7 as function of the mean value of the window position q. Note the different convergence rate 
for every a. All units are arbitrary. 



The values of the phase diffusion coefficient D are com- 
puted according to Eqs. (fI3]l, (fI5]l. and ([2D]). As for the 
solution of the adjoint operator, in these equations <j>- 
derivatives are calculated through Fourier differentiation, 
whereas y-dcrivatives are discretized by a fourth order 
centered finite differences, i.e. 



d q u (q) = 



1 



12A<7 



- u (q + 2Aq) + 8u Q (q + Aq) 
-8u (q-Aq)+u {q-2Aq) 



(31) 



These four additional stationary solutions, which are 
equally distributed around uo(<p;q) and separated by a 
factor Ag, are found by keeping constant the number 
of discretization points but changing the length of the 
system A, that is 2n/q. In this manner we are able to 
directly compare the values of these four functions and 
readily obtain the discrete representation of d q uo(4>;q). 
The integrals involved in the calculation of D are evalu- 
ated numerically by the extended Simpson's rule [l9j ). 

The functional relation between the phase diffusion co- 
efficient and the wavelength is found by computing D at 
discrete increasing system size (keeping constant the spa- 
tial discretization in cf> space). After the estimation of D 
for a given qi , the stationary solution for q2 < qi is com- 
puted by using the information provided by ito(gi)- The 
starting guess we used to determine 1*0(92) was a vector 
of N2 — 2-k I qiA<\) Fourier modes formed by the Fourier 
transform of uo(q\) for the first N\ = 2it /q\A<j) elements 
and zero padded for the others N2 — N\. Moreover, the 
amplitude of this initial guess has been adjusted accord- 
ing to the previously observed solutions. In fact, during 
the estimation of uq for decreasing values of q, we record 
the amplitude of these solutions and we are able to fore- 
cast the amplitude of the successive stationary solution 
(at least roughly). In this manner we are able to reduce 



the computational cost needed to estimate the new value 
of D. Finally, after the estimation of D for a set of in- 
creasing wavelength, the coarsening law follows by the in- 
version of relation (fig]), which leads to L(t) ~ y/\D(L)\t. 



IV. RESULTS 

We now report and discuss the results we obtained 
using the numerical methods described in previous Sec- 
tions. Let us start with the simplest and well known 
models: the GL model (|4aj) and its conserved version, the 
CH model (|5a|) . They admit analytical solutions based 
on the Jacobi elliptic function Sn(x;p) (24|, allowing us 
to test our numerical procedure with a controlled result. 
This family of solution is parametrized by the elliptic 
modulus p G [0, 1] (for more details see [E]) and take the 
form 



( AK 

uo(x) = A Sn I —j-x;p 



(32) 



The dependence of the complete elliptic integral of the 
first kind K and the amplitude A on the elliptic mod- 
ulus p is reported in (|B4| and (|B5[) . respectively. It 
has already been shown that for these two models the 
value of the phase-diffusion coefficient can be expressed 
as a function of A, the amplitude of uq and the integrals 
J = Jq(u' ) 2 dx for the GL model or / = J* Uq dx for its 
conserved counterpart 0, 0] 



D 



D 



A 2 (A - A 3 ) 
Jd A X : 

X 2 (A - A 3 ) 
Id A \ ■ 



GL equation 



CH equation 



(33) 
(34) 



As depicted in [B] we can compute all these functions 
and their derivatives with respect to A through complete 
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FIG. 7: (a) Absolute value of the phase diffusion coefficient D of the cKS equation ([6| in the case of r = for different values of 
A. We have used a fixed mesh discretization A</> = 0.025 and a variable q discretization Aq — 1.0 x 10~ 3 , 2.5 x 10~ 4 , 1.0 x 10~ 4 
ensuing a ratio Aq/q < 10~ 2 for the different values of A. (b) Dynamical evolution of equation ((6| for r = with simulation 
parameters N = 2048, Acj> = 0.1, and At = 10 -5 . The values of A are averaged over 50 different random initial conditions 
whose elements are sampled from a Gaussian distribution with zero mean and unit variance. The black solid line is a guide to 
eyes with slope 1/2. All units are arbitrary. 



elliptic integrals of first or second kind and Jacobi ellip- 
tic functions. Finally, the integrals / and J are easily 
evaluated by means of adaptive quadrature so that the 
values of D arc obtained without solving any differential 
equation. 

In the absence of noise, the Ginzburg-Landau and the 
Cahn-Hilliard model display a logarithmic coarsening, 
L(t) ~ C _1 lni, which corresponds, in the formalism of 
phase diffusion equation, to a diffusion coefficient which 
decreases exponentially with the wavelength of the sta- 
tionary solution, D(X) ~ e~ CA . Figure [4jt shows such 
exponential behavior on a lin-log scale while Figure 3)d 
shows the convergence of the constant C to a value of 
order 0.75 for both models. Note that in Figure 0Ji the 
comparison between the values of D estimated numer- 
ically (symbols) and those obtained through Eqs. ([55)1 
and (lines) confirms the validity of our method in 
a well-controlled situation and its applicability to other 
non-trivial cases. 

Next, we consider the nonconserved and conserved 
q— models, given by Eqs. (|4b"f and (|5b| . respectively. Ac- 
cording to the analytic treatment of the phase diffusion 
equation [l| and to numerics (taking care of dangerous 
caveats [25(), we expect that the coarsening exponent for 
the nonconserved models is n = 1/2 for 1 < a < 2 and 
n = a/(3a — 2) for a > 2. Case a = 2 has a loga- 
rithmic correction, L(t) ~ ^/t/lnt. Conserved models, 
instead show a constant coarsening exponent, n = 1/4, 
for all a. Figure[SJi considers cases a = 1,2 for non- 
conserved models, therefore we expect a constant D for 
a = 1 and D rj 1/lng for a = 2. This is exactly what 
our numerics shows. In Fig. [SJd we consider a — 1,2,3 
for the conserved models and, as expected, \D\ ~ q 2 for 
small q. The conserved case is analyzed in more details 



in Figs. |B^,b. On the left we plot the exponent 7 defined 
by \D\ ~ q 1 , and on the right we specialize on its con- 
vergence to the asymptotic value 700 = 2. According to 
the results displayed in Fig. |BJd, for a larger exponent a 
we observe a slower convergence to the asymptotic 700. 
This effect is more pronounced for the case a = 3 where 
a crossover region between two power laws has been al- 
ready reported in ;25|. Taking into account that this 
region is of the order A ~ 10 3 , i.e. q ~ 0.006, and for this 
wavelength our algorithm estimates a \D\ smaller than 
10 -5 , it is difficult to assess a clear asymptotic value for 
7 and the relative error \y— 7oo|/7oc saturates to a value 
around 5%. 

In Figs. [7^,, b we consider the cKS model ([5]). Some nu- 
merical and approximated analytical estimations of the 
coarsening exponents gave n — 1/2. i.e. a constant phase 
diffusion coefficient. This is shown on the left, while on 
the right we shows a clear evaluation of the coarsening 
exponent through the direct integration of the dynami- 
cal equation. Our results are much cleaner of available 
results in the literature. 

Finally we compare the computational cost used to es- 
timate the coarsening exponent by means of the phase 
diffusion method with the cost of the standard method 
used in the literature, that is, typically, the direct inte- 
gration of the time dependent PDE. We consider three 
cases of conserved models that are usually stiffer that 
their non-conserved counterparts due to high order spa- 
tial derivatives. Lines in Fig. [S] are found after the inver- 
sion of the relation A 2 /|Z?| ~ t that leads to the coarsen- 
ing law L(t) ~ t n . In the same figure, symbols are data 
obtained from the numerical integration of time depen- 
dent PDEs averaged over several realizations of differ- 
ent initial conditions (see captions for more details). As 
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FIG. 8: Comparison between the phase diffusion method and the numerical integration of the time dependent PDE for three 
different models, (a) Case of logarithmic coarsening, the Cahn-Hilliard equation. The black line is the estimation of A as 
function of t from the data displayed in Fig. |4ji. The red circles are the numerical integration of Eq. (|5a[) by means of our 
IFABM4 code with parameters L = 102.4, Ax = 0.1, At = 0.05, and averaged over 50 different initial conditions (generated 
from a normal distribution with zero mean and unit variance), (b) Case of power-law coarsening. As for the left panel, lines 
stand for the estimation of A as function of t from the values of D(X) while the symbols are obtained from the numerical 
integration of time dependent equations. The black solid line (from data in Fig. [SJa) and the blue circles are data for the 
conserved a-model Eq. (|5b|) with a = 2 and parameters L — 256, Ax — 0.25, At = 0.05, and averaged over 50 different initial 
conditions. The red dashed line and the green rhombi correspond to the conserved Kuramoto-Sivashinsky equation ©. These 
data are computed (line) and reported (rhombi) from the left and the right panel of Fig. [7] respectively. All units are arbitrary. 



shown in the figure, the phase diffusion method allows 
one to extract the coarsening law for times that arc orders 
of magnitude larger compared to the standard method. 
Moreover, the times listed in Tab. [TT| demonstrate the 
remarkable performance of our implementation of the 
phase diffusion method. Note that in the three cases we 
have not employed the same strategy to find the station- 
ary solutions of the different PDEs. As already explained 
in the previous sections, for the Cahn-Hilliard model we 
have only used our IFABM4 code with variable time 
stepping, for the conserved Kuramoto-Sivashinsky equa- 
tion we have solved directly Eq. ([23]) with the Newton- 
Raphson method, and for the conserved a-model we com- 
bined these two methods. We have also verified that the 
solution of the adjoint problem (flU]) for the cKS model is 
much less costly (roughly two orders of magnitude) than 
the determination of uq. 



V. CONCLUSIONS 

In this work we have demonstrated that the phase dif- 
fusion method can be a reliable and fast approach to find 
the coarsening law of nonlinear ID partial differential 
equations. Its numerical implementation is straightfor- 
ward and extends the applicability of the method beyond 
the cases studied in [l[: in fact, all these cases allowed an 
analytical evaluation of L(t), which is not always possible 
(see the cKS equation). 

Our algorithm permits to evaluate the coarsening ex- 



ponent with a negligible computational cost compared to 
standard methods, such as the direct numerical solution 
of the time dependent PDE (sec Table 2). Furthermore, 
it does not suffer from spurious effects associated with 
the finite size of the system or from errors arising in the 
time discretization of the PDE, and the result does not 
need to be averaged over several initial conditions. 

The main limitation of the method is the finite accu- 
racy in the estimation of D, which should be improved 
for larger times, when coarsening is slow (logarithmically 
slow or following a power law with a small exponent n). 
However, in such cases the direct time integration would 
require to attain extremely large times. Our discussion 
and in particular Table 2 shows that the phase diffusion 
method is computationally much faster than time inte- 
gration. 

A separate comment should be made for PDE whose 
steady solutions are determined by an ODE of higher 
order than two. In these cases |26| . which include the 
Swift-Hohcnbcrg equation, we expect more branches of 
stationary solutions, Xi(A), which makes our analysis 
more complicated. 

Finally, in this paper we confined our discussion to one 
dimensional PDE, but we have recently shown [27} that 
the phase diffusion method can be successfully applied to 
coarsening processes in 2D as well. It will be interesting 
to extend our numerical approach to the bidimensional 
case, specially for studying those equations which can 
hardly be attacked analytically. 
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Phase diffusion 
(long time) 


Integration 
(short time) 


Integration 
(long time) 


Ratio (long time) 
Integration/Phase diffusion 


CH 

ca-model (a=2) 
cKS 


8.32 x 10 2 
1.66 x 10 3 
6.07 x 10 2 


8.54 x 10 4 
3.95 x 10 4 
1.46 x 10 6 


~ 4.21 x 10 9 
~ 7.41 x 10 9 
~ 1.47 x 10 9 


~ 5 x 10 6 
~ 4 x 10 6 
~ 2 x 10 6 



TABLE II: Actual computational time (in seconds) required to measure the coarsening exponent n with the phase diffusion 
method (second column) and with the numerical integration of the time dependent PDE (third column). The fourth column 
is an estimation of the computational cost needed by the standard method to reach the same time span of the phase diffusion 
method for the data reported in Fig. [5] The fifth column is the ratio between the estimated time for the standard method 
and the cost for the phase diffusion method. The numerical experiments were carried out in a machine equipped with a single 
processor Intel® Core™ 2 Duo with a clock frequency of 2.4 Ghz. 
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Appendix A: Pseudospectral algorithm 

The solution of the linear part of Eq. ([21]) suggests the following change of variable 
z k {t) = e- Ukt u k (t), 

and the temporal evolution of this new variable depends only on the nonlinear operator 



(Al) 



d t z k (t) 



N [u{t)[ 



h ■ 



(A2) 



This is our starting point to discretize Eq. (|A2[) in the time domain by using one of the standard algorithms available 
in the literature. In order to enhance the stability properties of our method we have chosen a fourth order predictor- 
corrector method, the Adams-Bashforth-Moulton (IFABM4) scheme. The predictor step applied to Eq. (|A2|) is 



z k (t + At) = z k {t) 



At 
24 



55e 



-ui k t 



N[u(t)] k -59 



,-ui h (t-At) 



N[u(t-At)] k 



+ 37 e -«»(t-2Ai) N _ 2At)} k - g e -«»(*-3At) N ^ _ 3A ^ 



(A3) 



that leads to a prediction for u k at the next time step 



u k (t + At) = e UkAt 



«*(*) + ( 55N o - 59Nr + 37N 2 - 9N 3 ) 



(A4) 



where the functions N„ = e nuikAt N [uk(t — nAt)] k are computed during the time stepping of the algorithm by succes- 
sive multiplications of the factor exp(wfcAi). As before, the corrector step is applied to the variable Zk(t + At) and is 
computed from the value of the predictor and the value of Uk at three previous time steps 



z k {t + At) = z k (t) + 



At 

24 



9e -^( t +At) N [u(t + At)} k + 19e-" fct N [u(t)] k 
- 5e-^ ( *- At) N [u(t - At)] k + e -^(t-2At) N ^ t _ 2At)] k 



(A5) 



so that the update of u k reads 



u k (t + At) = - At N [u(t + At)], 



Ufc(t) + ^ (19N - 5Ni + N 2 ) 



(A6) 
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Equations (|A4[) and (|A6[) depend on four previous values of the function Uk, we have to compute the four initial time 
steps of the evolution of (f2~Tj) by an integrator of the same accuracy, for example a classical fourth order Runge-Kutta 
method (RK4). For the variable Zk the RK4 method reads 



At 

z k (t + At) = z k (t) + — (/i + 2/ 2 + 2/ 3 + U) , 
o 



(A7) 



where the different RK4 steps are obtained through a sequence of new variables. The first point in the RK4 evolution 
is easily found by evaluating Eq. (|A2|) at time t with the variable Uk(t) 



ft = e~" kt N 



T- 1 [e»« t z k (t)] i 



e^ fct N . 



Then we introduce the new variable u\k an d we find the second RK4 point ji 



11 1 = J' 1 



U k {t) 



At 



■N 



hence f2 = cxp[— ujk(t + At/2)] N [ui] k . The same procedure is iterated, so that 



u 2 = J' 1 



a u k At/2 



so that fz = exp[— + At/ 2)] N [112} k , and 
u 3 = F~ 1 [e^ Af (u k (t) + AtN[u 2 ] k )]^ 
hence fn = exp[— ujk(t + At)] N [1*3] fc . Finally, the update for Uk reads 



u k {t + At) 



M fe (t +ir N o 
o 



At 
~6~ 



2e u k At/2 (N[ui]fe + N[ua]fc) + N[M3]fe 



(A8) 



(A9) 



(A10) 



(All) 



(A12) 



This last equation is used to warm up the predictor-corrector method with the three initial time steps. 



Appendix B: Analytical solutions of the Ginzburg-Landau equation 

The stationary states of the Ginzburg-Landau equation can be expressed analytically by means of the Jacobi elliptic 
function sine-amplitude Sn(x;p) [23 |. This family of solutions is parametrized by the elliptic modulus p <E [0, 1] and 
take the form 



u (x) 



2p 2 



Sn 



P 2 + 1 \^/p T TT 



and its derivative reads 
u' (x) = - Cn 



P 2 + 1 V v^ + T 



p Dn 



where Cn and Dn are the other two Jacobi elliptic functions [28[. The periodicity of these functions is 

A = 4A»vV + l, 
where the complete elliptic integral of the first kind is defined according to 

dt 



^{l-t 2 )[l-pH 2 ] 

Taking into account that the function Sn G [—1,1], the amplitude of uq is a function of the elliptic modulus 



A = 



I 2p 2 

p 2 + r 



(Bi) 



(B2) 



(B3) 



(B4) 



(B5) 
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so that we are able to express the periodicity of the stationary solution as a function of its amplitude 

A 



X(A) = 4 



K 



2 -A 2 \V2~A5, 
The variation of A with respect to the amplitude of Uq is 



A \ 2 - A 2 



1/2 



1 



E 



A 



I -A 2 \V2~A^ 



K 



A 



V2~A^ 



where the complete elliptic integral of the second kind is 



E(p) = / dt 



2+2 



1-pH 
1-t 2 



(B6) 



(B7) 



(B8) 
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